NMDS Class Demo

Package installation and declaration (Chunk 1)

# vegan contains all the main nmds functions, as well as the data set used
if(!require(vegan)){install.packages("vegan")}

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.5/permute_0.9-8.zip'
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.5/vegan_2.7-2.zip'

package ‘permute’ successfully unpacked and MD5 sums checked
package ‘vegan’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\jedrake\AppData\Local\Temp\RtmpuuwHTe\downloaded_packages

Basic Example Visualization (Chunk 2)

NMDS on Basic Example (Chunk 3)

set.seed(1) #ensures reproducibility
vegdist(fourD) # uses bray-curtis distance by default, this is a step in metaMDS and doesn't need to be done here, but good for visualization

fakeNMDS <- metaMDS(comm = fourD, distance = "bray", trymax = 40)

ordiplot(fakeNMDS,type="n")
orditorp(fakeNMDS,display="species",col="red",air=1)
orditorp(fakeNMDS,display="sites",air=1)

Load and visualize data (Chunk 4)

# BCI is a data set in the vegan package that has data on tree counts from plots on Barro Colorado Island in Panama
data(BCI)
data(BCI.env)

# Combine and view environmental conditions data frame with species count data frame
head(BCI)
head(BCI.env)

# For viewing the entire dataframe
#View(BCI)
#View(BCI.env)

#Run NMDS (Chunk 5)

set.seed(67) # NMDS is iterative and random, so seed is set to ensure reproducible results (ensures the demo looks the same for everybody)
BCI_NMDS <- metaMDS(comm = BCI, distance = "bray", try = 40, trymax = 100) # the one command that does everything
print(BCI_NMDS)

Stressplot (Chunk 6)

stressplot(BCI_NMDS) # create stress plot visualization (as close as we come to checking conditions)

Basic Plotting options, works better for smaller data sets (Chunk 7)

# plots ordinations plainly, red dots are species representations of distance while grey dots are plots
# We can think of this as each species representing a vector (from the origin of the plot to the species dot) that pulls the plot dot towards it
plot(BCI_NMDS)

# basic plot option displaying most information possible, but ugly
ordiplot(BCI_NMDS,type="n") #Creates empty plot
orditorp(BCI_NMDS,display="species",col="red",air=1) #Writes in ordination vectors in red for "species." In this instance species is actually species (since this is the most common use of nmds), but in using for other situations, this could just be whatever the "dependent" variable is
orditorp(BCI_NMDS,display="sites",air=1) #plots points for each of the sites 

Clean Ordination Plots (Chunk 8)

# create array of habitat type for sorting plots 
hab <- as.array(BCI.env$Habitat)

#create array matching unique colors to each habitat type
colors = rep("",length(hab))
for (i in 1:length(hab))
{
  if (hab[i] == "OldHigh")
  {
    colors[i] <- "red"
  }
  else if (hab[i] == "OldSlope")
  {
    colors[i] <- "orange"
  }
  else if (hab[i] == "OldLow")
  {
    colors[i] <- "purple"
  }
  else if (hab[i] == "Young")
  {
    colors[i] <- "green"
  }
  else if (hab[i] == "Swamp")
  {
    colors[i] <- "blue"
  }
}

# create new plot scaled to necessary size for visibility, these numbers for xlim and ylim are guess and checked, plot must stay in correct aspect ratio to properly display distance
ordiplot(BCI_NMDS, type = "n", ylim = c(-.4, .4), xlim = c(-.8,.8)) |> 
  points("sites", col = colors, pch = 16) #plot points on this created plot for each site, with the designated color (pipe format)
ordihull(BCI_NMDS,groups=hab,draw="polygon",col=c("red","purple","orange")) # visually connect points with polygons

# Spider plot option
ordiplot(BCI_NMDS, type = "n", ylim = c(-.4, .4), xlim = c(-.8,.8))
ordispider(BCI_NMDS,groups=hab,col=c("red","purple","orange","blue","green"))

#NMDS Without Transformations (Chunk 9)

set.seed(67)
BCI_NMDS_noadjust <- metaMDS(comm = BCI, autotransform = FALSE, distance = "bray", try = 40, trymax = 100) #this line is the only one changed beyond an object name 
print(BCI_NMDS_noadjust)
stressplot(BCI_NMDS_noadjust)
ordiplot(BCI_NMDS_noadjust, type = "n", ylim = c(-.4, .4), xlim = c(-.8,.8)) |> 
  points("sites", col = colors, pch = 16) 
ordihull(BCI_NMDS_noadjust,groups=hab,draw="polygon",col=c("red","purple","orange")) 
ordiplot(BCI_NMDS_noadjust, type = "n", ylim = c(-.4, .4), xlim = c(-.8,.8))
ordispider(BCI_NMDS_noadjust,groups=hab,col=c("red","purple","orange","blue","green"))